{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Create climatology and persistence forecasts\n",
    "\n",
    "In this note book we will create the most basic baselines: persistence and climatology forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from src.score import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "sns.set_style('darkgrid')\n",
    "sns.set_context('notebook')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Load data\n",
    "\n",
    "First, we need to specify the directories and load the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "res = '5.625'\n",
    "DATADIR = f'/data/weather-benchmark/{res}deg/'\n",
    "PREDDIR = '/data/weather-benchmark/predictions/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# Load the entire dataset\n",
    "z500 = xr.open_mfdataset(f'{DATADIR}geopotential_500/*.nc', combine='by_coords').z\n",
    "t850 = xr.open_mfdataset(f'{DATADIR}temperature_850/*.nc', combine='by_coords').t.drop('level')\n",
    "data = xr.merge([z500, t850])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# Load the validation subset of the data: 2017 and 2018\n",
    "z500_valid = load_test_data(f'{DATADIR}geopotential_500', 'z')\n",
    "t850_valid = load_test_data(f'{DATADIR}temperature_850', 't').drop('level')\n",
    "valid_data = xr.merge([z500_valid, t850_valid])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Persistence\n",
    "\n",
    "Persistence simply means: Tomorrow's weather is today's weather."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "def create_persistence_forecast(ds, lead_time_h):\n",
    "    assert lead_time_h > 0, 'Lead time must be greater than 0'\n",
    "    ds_fc = ds.isel(time=slice(0, -lead_time_h))\n",
    "    return ds_fc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "lead_times = xr.DataArray(\n",
    "    np.arange(6, 126, 6), dims=['lead_time'], coords={'lead_time': np.arange(6, 126, 6)}, name='lead_time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "persistence = []\n",
    "for l in lead_times:\n",
    "    persistence.append(create_persistence_forecast(valid_data, int(l)))\n",
    "persistence = xr.concat(persistence, dim=lead_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:    (lat: 32, lead_time: 20, lon: 64, time: 17514)\n",
       "Coordinates:\n",
       "  * lead_time  (lead_time) int64 6 12 18 24 30 36 42 ... 90 96 102 108 114 120\n",
       "  * time       (time) datetime64[ns] 2017-01-01 ... 2018-12-31T17:00:00\n",
       "  * lat        (lat) float64 -87.19 -81.56 -75.94 -70.31 ... 75.94 81.56 87.19\n",
       "  * lon        (lon) float64 0.0 5.625 11.25 16.88 ... 337.5 343.1 348.8 354.4\n",
       "Data variables:\n",
       "    z          (lead_time, time, lat, lon) float32 dask.array<chunksize=(1, 8760, 32, 64), meta=np.ndarray>\n",
       "    t          (lead_time, time, lat, lon) float32 dask.array<chunksize=(1, 8760, 32, 64), meta=np.ndarray>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "persistence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "The forecast files have dimensions `[init_time, lead_time, lat, lon]`. Let's now save these files so we can evaluate them later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# Save the predictions\n",
    "persistence.to_netcdf(f'{PREDDIR}persistence_{res}.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Climatology\n",
    "\n",
    "First let's create a single climatology from the entire training dataset (meaning everything before 2017)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "def create_climatology_forecast(ds_train):\n",
    "    return ds_train.mean('time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "train_data = data.sel(time=slice(None, '2016'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "climatology = create_climatology_forecast(train_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (lat: 32, lon: 64)\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -87.19 -81.56 -75.94 -70.31 ... 75.94 81.56 87.19\n",
       "  * lon      (lon) float64 0.0 5.625 11.25 16.88 ... 337.5 343.1 348.8 354.4\n",
       "Data variables:\n",
       "    z        (lat, lon) float32 dask.array<chunksize=(32, 64), meta=np.ndarray>\n",
       "    t        (lat, lon) float32 dask.array<chunksize=(32, 64), meta=np.ndarray>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "climatology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEMCAYAAABkwamIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df3QU9b0//ufuJJuY4BKCJmzADyk5xa6XVs+J/eK9F7TGckGNRG5PDzZg9R6sXpEclZurASlgECVqrxUF6Q+0F00tl+tFMAWFXqhotXxsb9uLxFq/mCgmGyIJMUjAZGfenz8CW/Y9k5mdyc7ukH0+OHsOs+/58dr3zO4rM/va9/iEEAJEREQe4E93AERERGcwKRERkWcwKRERkWcwKRERkWcwKRERkWcwKRERkWcwKRERJZnQetIdwjnLlym/U1rwyC/QeezzdIeRmXzpDuAclhHvTu8oGjMKG5fclJR1RbvmAlrH0DP4xyFr7OakbGskyUp3AKnSeexzRLp60x1GZmJSco5J6Zylqu2A2jb0DIqWOR/ANrBPiIhcIE7/G4qPf3EYYlIiInKBBgEBbch2JiVjTEpERC6ICg2aGDop+U3aMhmTEhGRC1QIaCZnQ2aX9jJZxiQlNccHNTfBb9wdHCs+q2XsrjOB+S23mSZC7marbpfadcsnaxmT5R2x6P+E9o80j24Zq+kElvHJBbZuHDcW/Sl8BjNY7UOLadv73GidEjUneVU5mkVSYhWLMU/9Tmnv3r248cYbUVVVhRtuuAG7du0CALS0tGDu3LmYOXMm5s6di9bW1vQGSkRkQRMCqslDy4xf49jmmTMlIQTuu+8+NDY2YvLkyfjzn/+M73znO/jmN7+JFStWoLq6GlVVVdi2bRuWL1+OTZs2pTtkIqIhaacfQ+EvJYx56kzJ7/fj+PHjAIDjx4+jqKgIx44dQ3NzMyorKwEAlZWVaG5uRnd3dzpDJSIypUJYPkjPM2dKPp8PP/zhD7Fw4ULk5eXhxIkT+NGPfoRIJILi4mIoigIAUBQFRUVFiEQiKCwsTHPURETGomLwMRRevTPmmTOlaDSKH/3oR1i/fj327t2LZ555Bvfeey/6+vrSHRoRkW0qfJYP0vPMmdJ7772Hzs5OlJeXAwDKy8tx3nnnIScnB0eOHIGqqlAUBaqqorOzE6FQKM0RExENTRODD7N20vNMUho3bhw6Ojrw4YcfYtKkSTh06BCOHj2KiRMnIhwOo6mpCVVVVWhqakI4HLZ96W4gH+jvdxabG6W9ttuN5rH67V0i63SBVWmvrt1v0W7wB6XupdjdpoM/Uq1KkJNynGjm7Ub7XH7Op8rtPmlabjcI3Oq16Po3/gmhwHQaADSLeYR8HSeRkvBhnnxE84a3/Nk0i7MhP8+UDHkmKV144YVYuXIl7r77bvhOH+CPPPIICgoKsHLlStTV1WH9+vUIBoNoaGhIc7REROasLtExKRnzTFICgNmzZ2P27Nm658vKyrBly5Y0RERE5ExU+DGgO937K59JWybzVFIiIhopVPihmtSSmbVlMiYlIiIXDBY6DH2JjoUOxpiUiIhcYFXooPE7JUNMSkRELlDhh2ryvREv3xnLmKQUzQcGoqcnbJYPJ1Tqa1HKa1XObTVtuA2rZZyMNJ6M0c5t9q9l6a/Re3e4I0gn8EeqZWm7RdyOWB0nUrk3APij0jzStF8uEde16wO3eyzJfZFQSbj06aNlmy8jl5Ab9vcwfwoQzTdvt0ODH5pJ4jFry2QZk5SIiFJpQPjRb5SNT/Oz+s4QkxIRkQs0+Ey/N+J3SsaYlIiIXKBZlITz8p0xJiUiIheowqLQgZfvDDEpERG5gIUOzjApERG5QBOAyh/P2pYxSWkgqGHgdJ2r8Mt1txYL68qmDUpo7Y76rcm1q9L8RqNBS6W7liXhVmXqRttNwqjVVmyPIu5gNOhklITr1mlR9qxlCdN2AIAizWN35HGD8m25JNw/4JOm5XZpG9K00TotfzpgUS6vK+cGIOSScPnTSO5vi/432o6Q94l8ciLtj4HzHRzQQxgQWRiQX6TUTnrsFSIiF7DQwRkmJSIiF6jCZ3r5zqwtkzEpERG5YPB3SmZnSkxKRjyVlL744gs8/PDDePvtt5GTk4PLLrsMq1atQktLC+rq6tDT04OCggI0NDSgtLQ03eESEQ1JsygJ11gSbshTSemxxx5DTk4OXnvtNfh8Phw9ehQAsGLFClRXV6Oqqgrbtm3D8uXLsWnTpjRHS0Q0tAGhYMBkmCGztkzmmaR04sQJvPzyy3j99ddjt0O/4IIL0NXVhebmZjz33HMAgMrKSqxatQrd3d0oLCxMfAOjB4CsfgCAX66+syLNLoyuBVvMo9uiVUWfQYhCqtiTq/GsKvrkSjoA8FlVAcoVf/LgnkaDecrzWAzumdCAtxZsX55PZEBWy+ovqbIry7wdAERAWiY7vnN8uuo8+cDSB67K+2ggPlBfv1yNJ033G1T0yRV6dgf7TWBQXetBXOM3IlfnyZV1ACCyLfpXXkb6LBD5UtnhMAzeuiK5l+8qKioQCASQk5MDAKitrcX06dOHvMoEwPRKk9M2N3kmKR0+fBgFBQV4+umnsX//fuTn5+Puu+9Gbm4uiouLoSiDR6yiKCgqKkIkErGXlIiIUkiDz/wmfw6/U1q7di0mT54c99xQV5kA8ytNTtvc5JmLmtFoFIcPH8Yll1yC//qv/0JtbS1qamrQ19eX7tCIiGw7czt0s0cynLnKdPfdd8ddZQIQu9JUWVkJYPBKU3NzM7q7ux23uc0zZ0olJSXIysqKdcKll16KMWPGIDc3F0eOHIGqqlAUBaqqorOzE6FQKM0RExENTQi/aTGDON0WiUSgqvHXvIPBIILBoOFytbW1EEKgvLwcixcvRnt7u+FVpssvvxyRSGTIK01CCEdtbl+h8syZUmFhIaZOnYrf/OY3AAavZ3Z1daG0tBThcBhNTU0AgKamJoTDYV66IyJPU0/fDt3sAQDz5s3DNddcE/f493//d8N1NjY2Yvv27XjppZcghEB9ff2QV5k+//zzVL7cpPHMmRIAPPjgg1i6dCkaGhqQlZWFRx99FMFgECtXrkRdXR3Wr1+PYDCIhoaGdIdKRGQqKvymFXbR02dKjY2NhmdKRs5cIQoEAqiursadd96JpUuXGl5lamlpQUlJyZBXmoQQjtrc5qmkdNFFF+H555/XPV9WVoYtW7akISIiImc0i8t3Z9oS/aDv6+uDqqo4//zzIYTAjh07YleNzlxlmjZtWuwq08SJExEMBmNXmqqqqnRXmpy2ucknhMiIsWqrXn8UkZM9Cc0r94j8Azijihp5GcOycRuMlpe3a7VNy7J06MvM1Wj8X3Zqf/xrF7pyY/2bTi459kXl0vX4SauydEO6QVvNF0pkkFfd54euJFwecFXaplxunG0QU1b8i8/KiS9Bzs6O/4s5Oyt+2mfwOqOqIk3HBx7tj//b02qfAvqycnkgWLsD8Rp+Nsvl73L5ttR/fqlvsqRpAMiS+itL0UynFemFhM4rwPZv3GcQrH1Pvn8bPhvoHLJ9dHYR7r74pwmv7/Dhw6ipqYGqqtA0DWVlZVi2bBmKiopw+PBhLF26FD09PcjKysI999yDq666CgBw6NAh1NXVobe3N3aladKkScNqcxOTkgEmpbOmmZTOamdSirUzKVn64fvfQ49JUirILsI9F/8kKdsaSTx1+Y6IaKQYHJDV7M6zHPvOCJMSEZELNGHx41kmJUNMSkRELohajH0X5dh3hpiUiIhcoMFvcesKz/xM1FOYlIiIXKDC4iZ/vJ+SoYxJSv8neAx5OV0A9Pcxka/tRqX2qOY3bTdahyzLbrmSA36pMsvqdQFAv1S5dSqaHTfd1x8/fUqaHjilP4S0gfh16ir25K6wqr4z6lr5ObkyTpr26SrrDCrjpP6TR+z2S5VzilxJp8RXfgUMqsNys6LSdPxw3DlSe8AvVZP5rY8jeZ/2Sfv05EAgbvpUVL8P+6V9qErvAblqU0dqVgziDmTHv9bzdNP9cdOjpOlcRT+id5bUX/J7wsoFgeS9T4Uw/1zIjLpn+zImKRERpVKiP56leExKREQuGBxmaOjEY3TlgpiUiIhcwTMlZ5iUiIhcoMFneiM/pzf5G+mYlIiIXKAK81EbVBY6GGJSIiJyQaI3+aN4GZOUyvKOYmz24OCIp7T4EtkvtCzT6agWXx6byGl3lk8aWFOqg/brpq3/bFKkeeRtKNJIp/Ltlo1+Qf6FGv9aT0p9c3wgN266dyAnvr0/fhoATg5IZeNR8wFDZT6pe436RpEH1pRKjrP95uXaRiX68joCcom3VG4slyTnKvHl3ecp8SXMAJDjlwZglbbph/lxke3Tl5nL5DLkL0T8/vg8Gl8SflKNnwaAfumYl98DMjlOuTQ736Av8rO+iJvO80sl4Fmn4qbl127UF3ISkH8LJI+wMCC910cbDaLrEIcZciZjkhIRUSpF4TetsItyRAdDnuyVp59+GhdffDH+8pe/ABi8NfrcuXMxc+ZMzJ07F62trekNkIjIwpnLd0M9ePnOmOd65eDBg/jjH/+IkpKS2HMrVqxAdXU1XnvtNVRXV2P58uVpjJCIyJom/noJz/iR7gi9yVNJqb+/H/X19VixYgV8p79Y6OrqQnNzc+z+85WVlWhubkZ3d3c6QyUiMnWmJNzsQXqeSkpPPvkkZs+ejYsuuij2XCQSQXFxMRRl8AtKRVFQVFSESCSSrjCJiCyZnyWZF0FkMs8kpT/84Q84cOAAqqur0x0KEdGwCYuEJJiUDHmm+u6dd97Bhx9+iGuuuQYA0NHRgQULFmDJkiU4cuQIVFWFoihQVRWdnZ0IhUK21v+lnE9xQukAoC8LlUvETwm5ZFwqcU7g5lxyuWqOP75cON8fXw7rl4bGNjq1l+dRpLJm+dbLcrtR3PIy8mvv0+LLhT+P5pq2A8BJNX4dcim6/INCuVxeLic2Kv2VS6t1pdZy+bDUruiGIteX6cvrkEvu5fZcaR8bxZ3t049sfTZ9SbM0bfBBJv9UQCYvI9/Hx+iW3fLxZ3VsWZVr5/r0JeG6/oLU3xYj6xvFPQCp5FvEf8T1S8ei3N95SfxEjGp+3R0G5HbS80yv3H777XjzzTexZ88e7NmzB+PGjcPGjRtx3XXXIRwOo6mpCQDQ1NSEcDiMwsLCNEdMRDQ0fqfkjGfOlMysXLkSdXV1WL9+PYLBIBoaGtIdEhGRKWHxvREv3xnzbFLas2dP7P9lZWXYsmVLGqMhIrJHg/moDe7f9vPc5NmkRER0LuMwQ84wKRERuUDT/LrbyMvtpJcxSaksuxP9/sHfNslVaPKgjXIFmjx/v7DutlxffGVRQFedZF6FZVwRZT7YpLyMroLKoK5Fv4xFtZc0/qrROq1YVrFZ9N3gMnK1nfnFEHmbRpy8FrNtGFX4Wcdgvs8S+evaatDWROKyGiBYXoc87ZfCNKx2tIhB3mP6wVX1fSHf6dWqP+V9FlDi3/vDwfspOZMxSYmIKJV4+c4ZJiUiIhcIix/IsvrOGJMSEZELeKbkDJMSEZEbhMXZEEcJN8SkRETkAlX4oGpDJyWjIaOISYmIyBWsvnMmY5JSKEuFFiuVjS+ZVYU8GOqp+HbpPFtN4LQ74Is/4OwWGw8YnNvL27UqmZVvIia3G5HnsVv6CwDZuoFj5XZpG1JfZUu95TeIW15Gpt+n0j406F/N4rjQz+8+ozhlCsyPtWyfXAZtfRz4LY5Yq/5PhN19pGs32D/yPpTJr0p+HX7FfHk7WOjgTMYkJSKiVOLYd84wKRERuUCIwYdZO+kxKRERuYCX75xhUiIicoFqMfadWVsmY1IiInKBgMXlu5RFcm7xTKo+duwYvve972HmzJm44YYbsGjRInR3dwMAWlpaMHfuXMycORNz585Fa2treoMlIrIgxF8v4Rk/7K+zoqICs2bNQlVVFaqqqvDGG2/EtT/99NO4+OKL8Ze//CX2nNnnp9M2N3nmTMnn8+G2227D1KlTAQANDQ14/PHH8fDDD2PFihWorq5GVVUVtm3bhuXLl2PTpk221j/anwuI8xKaVxXDL/ZVfOb5Xt6GXO6qGRQcW5XQ6uZ3UE4sc1IaLC8jl3TLfWO1DSc0i1HDjfaxbh9I67AuEXf/b1+j8niZrszZ5v4w3q69fSQfv0b9nSW9lETeA3ES+ErG9mv351qvNFEW3ynB4XdKa9euxeTJk3XPHzx4EH/84x9RUlIS97zZ56fTNjd55kypoKAglpAA4LLLLkN7ezu6urrQ3NyMyspKAEBlZSWam5tjZ1FERF4kEngAQCQSwSeffBL36O3ttbWt/v5+1NfXY8WKFfCd9UeJ2een0za3eeZM6WyapuHFF19ERUUFIpEIiouLoSiD9zRSFAVFRUWIRCIoLCxMc6RERMaE5oMwGWboTNu8efPQ1tYW17Zo0SLU1NQYLldbWwshBMrLy7F48WIEg0E8+eSTmD17Ni666KK4ec0+P4UQjtrc/tz1ZFJatWoV8vLyMH/+fDQ3N6c7HCIi2xItCW9sbISqxo8kEQwGDZdpbGxEKBRCf38/Vq9ejfr6esybNw8HDhxAbW1t8oJPI88lpYaGBnz00UfYsGED/H4/QqEQjhw5AlVVoSgKVFVFZ2cnQqFQukMlIhpSoj+etfNZdmbeQCCA6upq3HnnnZg8eTI+/PBDXHPNNQCAjo4OLFiwAI888gjC4fCQn59CCEdtbvPMd0oA8MQTT+Ddd9/FunXrEAgEAABjx45FOBxGU1MTAKCpqQnhcJiX7ojI08wr7yyKIAz09fXh+PHjp9ctsGPHDoTDYdx+++148803sWfPHuzZswfjxo3Dxo0bMW3aNNPPT6dtbvMJ4Y3BLj744ANUVlaitLQUubmDFTATJkzAunXrcOjQIdTV1aG3txfBYBANDQ2YNGmSrfVHP70KUNusZwSr785mVe3l2eo7i8qthKrvdBVkrL5LVCLVd/plbFbfJcD2a1fGI+vC14e9XQC48pV1aDvx2ZDt4/NHY98NdyW8vsOHD6OmpgaqqkLTNJSVlWHZsmUoKiqKm6+iogIbNmyIVeiZfX46bXOTZ5KS27RPr044KaWC1Zs2kQ84N9609pe3X9aayIdgsjn5UNS3D7+/5cSWjNG2rfahV/eR1T5xkuTtJlxd3ynj4b9wr+3tGpm+3TopvTE78aSUKTz3nRIR0UiQaPUdxWNSIiJyw9k/RhqqnXSYlIiI3ODSiA4jHZMSEZEbeKbkCJMSEZFreDZkV8YkpQERBUTU0bLJqERyUl2XCnJFmf3SX/uvQxOq9UwuS0YlnRPJqLaTWb0WuTWRfWx3Hzmp8EvFe8Cqwk+Ve0dEkZOsjQvoO19uJ52MSUpERCklfObfG/E7JUNMSkRELkh0mCGKx6REROQGFjo4wqREROQGXr5zhEmJiMgFPjH4MGsnPSYlIiI3aL7Bh1k76WRMUhKn/znhhRJmIDVlzOkqlXab1QjfI5lchu7GQL7e+YnD8OLwJft1eKNbzikZk5SIiFKKhQ6OJPxLyZ07dxo+/+qrryYtGCKiEUMk8CCdhJPSAw88YPj88uXLkxaMmZaWFsydOxczZ87E3Llz0drampLtEhE5cqb6zuxBOpaX7w4fPgxg8Pa7Z/5/dtuZ25a7bcWKFaiurkZVVRW2bduG5cuXY9OmTSnZNhGRbRbVdzxTMmaZlGbMmAGfzwchBGbMmBHXdsEFF6Cmpsa14M7o6upCc3MznnvuOQBAZWUlVq1ahe7u7pTcM56IyLYM+07pu9/9Lm644QZ8+9vfjnv+9ttvx49//OOE12OZlP785z8DAObPn48XXnjBZpjJEYlEUFxcDEVRAACKoqCoqAiRSIRJiYg8KdN+p/SHP/wBXV1deO+99/DAAw/EPq9/97vf2VpPwt8ppSshESWD4vPFPYhcl2HfKWVnZ2Pz5s1oa2vDrbfeip6eHgCDX/3YkXBJeDQaxc9//nO88847OHbsWNyGGhsbbW3UrlAohCNHjkBVVSiKAlVV0dnZiVAo5Op2iYiGZYSdDVkZNWoUNmzYgH/7t3/Dt771Laxbtw4+m38EJnym9Mgjj2Dz5s24/PLLcfDgQfzDP/wDurq6cMUVV9gO3K6xY8ciHA6jqakJANDU1IRwOMxLd0TkXRlWEn7mRMXn8+Ff/uVfsHjxYvzTP/0T+vv7ba0n4aS0a9cu/OQnP8Ett9wCRVFwyy23YN26ddi/f7+9yB1auXIlXnjhBcycORMvvPACHnzwwZRsl4jICZ9m/RhJVq9eHTd9/fXXY+PGjfjnf/5nW+tJ+PLdqVOnYpfLcnNzcfLkSZSVlaG5udnWBp0qKyvDli1bUrItIqJhy7Dqu+uuu0733CWXXIJLLrnE1noSTkplZWU4cOAAvva1r2HKlCl46qmnMGrUKBQXF9vaIBFRJsi06rtkSfjy3dKlS5GVNZjD6urq0NzcjF//+tdYtWqVa8EREZ2zMqz6LllMz5TefvvtIZ+74447AAADAwMuhJV8vtP/vMLJaMbyyMz6dY6wi9QuYlm4OatjzXp5+/1r9Z5wsk77kriNDLt8lyymSWmo8e7O5vP58N///d9JC4iIaCTwweLyXcoiObeYJqU9e/akKg4iohHFqsJupFXfJQvvp0RE5AZevnOESYmIyA1MSo4wKRERuYAl4c4Mr8SGiIgoiTLmTEnx+QFfcnKwKob/DWUyylvlEtrhlvG6xQul6l7tG1k6+spJ37hRnp2Kkm/F6jMgSZ8RAHj5zqGMSUpERKnkExbVd0xKhpiUiIjcwDMlR5iUiIjcYFHowKRkjEmJiMgNPFNyhEmJiMgFbpSEV1RUIBAIICcnBwBQW1uLKVOm4L777sPHH3+MQCCAiRMnor6+PnYT1JaWFtTV1aGnpwcFBQVoaGhAaWnpsNrc5BN2b6DuggcffBBvv/02AoEA8vLy8MADD+CrX/0qAODkyZNYsmQJDh48CEVRcP/99+Pqq6+2vY1TndMg1DbDNjcq4dJRRaVa7MpEBoFVLeZRktBX8jbkuKxex+Ay8axqpuQBWJ3s82S8djcqzJwM7muXHLcXBrR1o6LSp4xHbtGbSVnXNx/diPae3iHbSwqC+NV9C2yts6KiAhs2bMDkyZNjz/X09OD999/H1KlTAQANDQ347LPP8PDDDwMAvvvd7+Jb3/oWqqqqsG3bNrz00kvYtGnTsNrc5Ik62SuvvBKvvPIKtm/fjjvuuAP33ntvrG3jxo3Iz8/H7t27sWHDBixbtgwnTpxIY7RERNbOnCmZPQAgEongk08+iXv09g6dzGQFBQWxhAQAl112Gdrb2wEAXV1daG5uRmVlJQCgsrISzc3N6O7udtzmNk9cvjv7zOeyyy5DR0cHNE2D3+/Hzp07sWbNGgBAaWkppkyZgn379uHaa69NV7hERIlJ4CR23rx5aGuLv4qzaNEi1NTUGM5fW1sLIQTKy8uxePFiBIPBWJumaXjxxRdRUVEBYDDhFRcXQ1EUAICiKCgqKkIkEoEQwlHbmcuCbvFEUjpbY2MjvvGNb8DvHzyJa29vx/jx42PtoVAIHR0d6QqPiCgxCRY6NDY2QlXVuKazE83ZGhsbEQqF0N/fj9WrV6O+vh6PP/54rH3VqlXIy8vD/Pnzhxl8+qQkKc2ZMyd2Oil76623Ytn4l7/8JV555RU0NjamIiwiItckWugQCoUSXueZeQOBAKqrq3HnnXfG2hoaGvDRRx9hw4YNsT/qQ6EQjhw5AlVVoSgKVFVFZ2cnQqEQhBCO2tyWkqS0detWy3l2796NJ554Aj/72c9wwQUXxJ4vKSlBW1tb7JQxEonEXT8lIvKkJJeE9/X1QVVVnH/++RBCYMeOHQiHwwCAJ554Au+++y5+/OMfIxAIxJYZO3YswuEwmpqaUFVVhaamJoTD4djnqdM2N3mi+m7v3r1YtWoVnnvuOUycODGu7amnnsKRI0fw0EMPobW1FdXV1di1axdGjRplaxusvmP1Xfz8rL6zg9V39s18aCPaj5lU340J4rVliVffHT58GDU1NVBVFZqmoaysDMuWLcNnn32GyspKlJaWIjc3FwAwYcIErFu3DgBw6NAh1NXVobe3F8FgEA0NDZg0adKw2tzkiaR0xRVXIDs7Oy4L/+xnP8OYMWPQ19eHuro6vPfee/D7/fjXf/1XfPOb37S9jU87/j+o6ieGbVYfWIl8GCXjg9aK/EEsb1OVNiHPPyD0r0OVXps8bUUx+EA0es5sGwMi/sOlH0rctJZA3FasYgIAv8UPRxSpR+V1yssn0jd+n3m7LkbT1kF2/xyy25eAweuQ21OQsxJ5X9pNW4oyAReO+7/OApLMXJVAUvq+vZLwTOCJQoff/va3Q7bl5eVh7dq1KYyGiGj4fKcfZu2k54mkREQ04nCYIUeYlIiIXOCDRfVdyiI5tzApERG5gWdKjjApERG5wKdZ3OQv/Tdk9iQmJSIiN/BMyZGMSUofR3PwhZqb0LxyuWu2TzVtHyRfIR5eqbVqUMwql0brS6vjd2e/iC+tHpBKrY2WUXXbsP9bEF3ptE8uXTffhlVMRvNocn+L4f+GRbH4U9YvHQfy/H6D4uyAdCz5pWXk9mxfNH4bBuuU47BL7jtAv0/k/gxIcVn1RSL7wygOO9sAjErV5WPRfJ/miBxcaDqHDbzJnyMZk5SIiFKKZ0qOMCkREbnAjZv8ZQImJSIiNwiYD6/BpGSISYmIyAU8U3KGSYmIyA38TskRJiUiIhf4hIDPZCBms7ZMljFJ6Z2Tk9A7kNjtLuQyXbkkXJ4G9CWyunVa/FlkVQ4LAP1SGfSAXPJtMf2Flq1b54BmvoxmURJuVPasHy1bM22XGY0KLrOK06iM3O42NJtl5fJxYyTbb142Lh9bcl8ZHXu6UnSbv8o0ep1WJdx2t5nINuyW9VuVdwMJlO1L06Ozi1BuudYE8UzJkYxJSkREqcTvlJxhUiIicoFPWAwzxKRkKPm3bhyG/fv3IxwO44UXXog9d/LkSdxzzz2YMWMGZs2ahb1796YxQiKiBIkEHqTjmTOlzz//HI8//tP/dCgAABI3SURBVDiuvPLKuOc3btyI/Px87N69G62trZg3bx527dqF/Pz8NEVKRGSNl++c8cyZ0po1a7BgwQKMGTMm7vmdO3fipptuAgCUlpZiypQp2LdvXzpCJCJKHM+UHPFEUnr99dfR29uLWbNm6dra29sxfvz42HQoFEJHR0cqwyMisu3MmZLZg/RScvluzpw5aG9vN2x79dVX8YMf/ADPPfecqzH86shX0HGqGIC+zFYuG/VJ7VlSGa88DQBZPvN5rMqgZUYly1Et/m+IqFQyK7cPqFK5t6YfJVyVlkmkVPpsikFfyM9l++PLmK36U94fdvvOiPy6EinB1y2ThDJzmdVr0x2rBvMPd5TwZJDjSqS/dcezZj4yubxOYdDfVvvAqr/H5Y7Bwi+bzpI4TcCnmWzPrC2DpSQpbd26dci23/3ud/j000/x7W9/GwBw7Ngx7N27Fz09PVi0aBFKSkrQ1taGwsJCAEAkEsHUqVNTETYRkXP8nZIjaS90uPzyy/H222/Hpuvq6jBlyhTMnz8fADBr1ixs3rwZX/3qV9Ha2ooDBw7gBz/4QbrCJSJKCEvCnfHEd0pmFixYgN7eXsyYMQN33HEH6uvrMWpUYiMzEBGlDQsdHEn7mZJszZo1cdN5eXlYu3ZtmqIhInKGJeHOeC4pERGNCEIMPszaSSdjklLH/38B2o4HEptZKuARfungMbroKc2jG0vS7p9FBlVEuuvTmjSPGj8tz28YgrwO3UalsKxeJwBkx29YZEnVjfK0IlXfZUnTcv8D8EnPyVVVcgWlzKhyS15CN48wb5c/Y4y2YbRfzdaRFE7WKYXpk7tC2GsXqkF/y8eeNC0v45PXYbBOeR26w8CiL3znB4Hp5vMkyqdZfKdkb9zcjJExSYmIKJV4+c4ZJiUiIldYXL5jpYMhJiUiIhfwTMkZJiUiIjfwx7OOMCkREbmAZ0rOMCkREblBFYMPs3bSyZikNPoDP/q6Ttcv2z0W5MpTo0pfqTRaV/lrd6xOgxh1Jd6W03JdrsE65Vl0JeBSma4iDbRpcARp2X7TeTSpMl/Llqal+YXRNpT4wKWxZ63732h/WJQPW5YXW81vFJj8awMn67SShJJwu8ez/ucL1vP4o1J71F47AMjjA/tUadriPTO6MHmD3PBMyZmMSUpERKnF6jsnmJSIiNxgdc8k5iRDTEpERG5g9Z0jTEpERC7wqYDPpJhB/r4rERUVFQgEAsjJyQEA1NbWYvr06WhpaUFdXR16enpQUFCAhoYGlJaWAoArbW7y/K0riIjORT4hLB9OrF27Ftu2bcO2bdswffrgQH0rVqxAdXU1XnvtNVRXV2P58uWx+d1ocxOTEhGRG1J0P6Wuri40NzejsrISAFBZWYnm5mZ0d3e70uY2z1y+e/7559HY2Ijs7GwoioKXX34ZAHDy5EksWbIEBw8ehKIouP/++3H11VfbXv/oQ6dwquMkAKPSXmmEaQej9+pHBZdKqW2WhBuXE0txypcGdOXD9o96IQ/3bFkSrv+7RsuW5gn4zduzfBbt+jg1uQRcF6e0gFWJswN2S8QdrcMrbPaf/Lp0P0+AdUm4PyqNBD9g3g7o3xP+AfP3jDwdLO7XrdO5xKrvIpEIVDX+Wl4wGEQwGDRcqra2FkIIlJeXY/HixYhEIiguLoaiDB70iqKgqKgIkUgEQoiktxUWFg6nUyx5Iint2rULr776Kv7zP/8To0aNwqeffhpr27hxI/Lz87F79260trZi3rx52LVrF/Lz89MYMRGRuUR/pzRv3jy0tbXFtS1atAg1NTW6ZRobGxEKhdDf34/Vq1ejvr4et956axKjTj9PJKVnn30Wd999d+w25xdeeGGsbefOnbG70ZaWlmLKlCnYt28frr322rTESkSUkARv8tfY2Gh4pmQkFAoBAAKBAKqrq3HnnXdiyZIlOHLkCFRVhaIoUFUVnZ2dCIVCEEIkvc1tnvhO6dChQ/jTn/6Em266Cf/4j/+I//iP/4i1tbe3Y/z48bHpUCiEjo6OdIRJRJQwnyosH8DgZ9qECRPiHkZJqa+vD8ePHwcACCGwY8cOhMNhjB07FuFwGE1NTQCApqYmhMNhFBYWutLmtpScKc2ZMwft7e2GbW+99RZUVUUkEsHPf/5zHDt2DN/5znfwpS99CV//+tdTER4RUfIl+XdKXV1dqKmpgaqq0DQNZWVlWLFiBQBg5cqVqKurw/r16xEMBtHQ0BBbzo02N6UkKW3dutW0vaSkBJWVlfD7/Rg7diz+7u/+Dv/7v/+Lr3/96ygpKUFbW1ssQ0ciEUydOjUVYRMROWZV9m23EOmiiy6KFYDJysrKsGXLlpS1uckTl+8qKyvxxhtvABg8Rf3973+Pr3zlKwCAWbNmYfPmzQCA1tZWHDhwIFabT0TkXeKv3ysZPTxbaplenih0uPXWW/H9738f119/PQCgqqoKf//3fw8AWLBgAerq6jBjxgz4/X7U19fHCiLsyO7oRaCtZ3BClWpRhTzt4GCRS6l9Ur73y+0WNbVGMagOfgJuFhNgHZfVtKJfp1w2DkUxb5fL5/1y3+k2YbCM1TrNlweSUyZuyWq/6+a3nkVfxm9vE7rlE9iuftRw+8ezviQ8/gld+bauXf/bDd9A/HtEXgbyOrT46YDt4fxNaDAcHT2unXQ8kZRyc3Px2GOPGbbl5eVh7dq1KY6IiGh4kn35LlN4IikREY04mgA0k9Mhgx8UE5MSEZE7ePnOESYlIiIX+GBx+Y6FDoaYlIiI3JDgiA4Uj0mJiMgNTEqOZE5SOvUFcPLU4P/lg8Hsy0ij+RMp65XLmu2WzBrFZBWHVHqti0EeOdvoSXkbVm8cgzJ1q97RjRitSevQlewbxOBknwyX1TYsRlg3nEcu05dL7K1K8g2WsSwRT2Sdckm9RRm/flT8+GnDy1iaeQm4rl0u7zYoCde9b+R5rI61vlz9Op1SxeDDrJ10MicpERGlktWN/HimZIhJiYjIDbx85wiTEhGRGwTMf4vEnGSISYmIyA08U3KESYmIyA1MSo5kZlLSDfhpUfFkVEUls1v9JR+QchWQ0frkZXTVdfLrkirrsgzK74ZbFWg0SKyu4kkeFDMavwlpWl5eDAxYb0MeVBdycwIfAPI6pMo4n27QV/P+98nVkIBBhaRUpZYlvSWtKioB3T7UxSkfB4kMDiyvUz+HPYlUUMr7yGqgZKPqu0S2m8z5zaiaeYyJxJ+BMjMpERG5TWjmfyxZ/CGVqZiUiIhcYXH5jpUOhpiUiIjcoMG8+o4nSoaYlIiI3MBCB0c8cTv0lpYW3HzzzaiqqsK1116Lp556KtZ28uRJ3HPPPZgxYwZmzZqFvXv3pjFSIqIEmd0K3SphZTBPnCk99thjmDlzJubPn48TJ06gsrISV111Fb72ta9h48aNyM/Px+7du9Ha2op58+Zh165dyM/PT3fYRERDU1Xj6tSz20nHE0nJ5/Ph+PHjAIBTp07B5/OhsLAQALBz506sWbMGAFBaWoopU6Zg3759uPbaa+1tJDoAnCktTmRAyrPJZbkJlNBalt3qynQtynaTIWrwJrBbduugJFzI88gl4XK7VWmw0TJyewLrsBa/Dd3ftVYl44mUhOsG1ZXK46X5jcvMrX7SYHNwYCes/uq3GvTYaB1OziTs3s1V3mdRg58fOMZCByc8cflu6dKl2LFjB6ZPn46KigosWLAAEyZMAAC0t7dj/PjxsXlDoRA6OjrSFSoRUWJ4+c6RlJwpzZkzB+3t7YZtb731FjZv3oyqqircdttt6OzsxM0334wpU6bg0ksvTUV4RETJx+o7R1KSlLZu3Wra/vzzz+NXv/oVAKCoqAhXXHEF3nnnHVx66aUoKSlBW1tb7HJeJBLB1KlTXY+ZiGhYhAbBH8/a5onLdxMmTMAbb7wBAPj888/x+9//Hl/+8pcBALNmzcLmzZsBAK2trThw4ACmT5+etliJiBJyZpghswfpeKLQ4ZFHHsFDDz2EZ599FtFoFNdddx2uuuoqAMCCBQtQV1eHGTNmwO/3o76+HqNGjUpzxEREFoRmXuDBMyVDnkhKU6ZMwS9+8QvDtry8PKxduzbFERERDRN/POuIJ5JS2lkdHPIo1m6U1KaL1Wu3LPVN4I0ll90GsuMmfZrFYWgQg+UeSKQEOdmMRvC2YvNYMvyOImr+Wn3yaObnqkR+JuHGTykcEpqAMDkOExq5PgMxKRERuYFnSo4wKRERuUETFiXhTEpGmJSIiFwgNNV05BGhcZghI0xKRERuEMLiJn88UzKSMUnpglDBXyeGXajgnS9Th2+4hQ6JbMJiJqvLGE7evOkot3VSUGD3WHTyRf6IKXRwfxNxnxPDNLZkjGkxw9iSMUnb1kjiE4LpmoiIvGGE/AlFREQjAZMSERF5BpMSERF5BpMSERF5BpMSERF5BpMSERF5BpMSERF5BpMSERF5BpMSERF5RkYMM9TS0oK6ujr09PSgoKAADQ0NKC0tTXkcFRUVCAQCyMnJAQDU1tZi+vTpKY2voaEBr732Gtra2vDKK69g8uTJAMz7yO34hoppqP5KRUzHjh3Dfffdh48//hiBQAATJ05EfX09CgsL09ZXZjGls68AYOHChfjkk0/g9/uRl5eH73//+wiHw2k9roaKKd19RRZEBrj55pvFyy+/LIQQ4uWXXxY333xzWuK4+uqrxfvvv697PpXxvfPOO6K9vV0Xi1kMbsc3VExD9VcqYjp27Jj47W9/G5tes2aNWLJkieW23YzLLKZ09pUQQvT29sb+v3v3bnHjjTdabtvtuIaKKd19ReZGfFI6evSoKC8vF9FoVAghRDQaFeXl5aKrqyvlsRi9GdIV39mxmMWQyvgSTUrp6LNXX31V3HLLLZ7pq7NjEsJbfbV161YxZ84cT/XVmZiE8FZfkd6Iv3wXiURQXFwMRVEAAIqioKioCJFIBIWFhSmPp7a2FkIIlJeXY/HixZ6IzywGIURa45P7KxgMprzPNE3Diy++iIqKCs/01dkxnZHuvnrggQfwm9/8BkII/PSnP/VEX8kxnZHuvqKhsdAhhRobG7F9+3a89NJLEEKgvr4+3SF5mlf6a9WqVcjLy8P8+fPTsn0jckxe6KvVq1fj17/+Ne699148+uijKd++EaOYvNBXNLQRn5RCoRCOHDkC9fQdIFVVRWdnJ0KhUFpiAYBAIIDq6mr8z//8jyfiM4shnfEZ9ZdVvMnW0NCAjz76CD/84Q/h9/s90VdyTIA3+uqMG2+8Efv378e4cePS3ldyTMeOHfNUX5HeiE9KY8eORTgcRlNTEwCgqakJ4XA45afjfX19OH78OABACIEdO3YgHA57Ij6zGNIV31D9ZRVvMj3xxBN49913sW7dOgQCActtpyIuo5jS3VcnTpxAJBKJTe/ZswejR49Oa18NFVNOTk7ajysylxE3+Tt06BDq6urQ29uLYDCIhoYGTJo0KaUxHD58GDU1NVBVFZqmoaysDMuWLUNRUVFK43vooYewa9cuHD16FGPGjEFBQQF++ctfmsbgdnxGMW3YsGHI/kpFTB988AEqKytRWlqK3NxcAMCECROwbt26tPXVUDHV1dWlta+OHj2KhQsX4uTJk/D7/Rg9ejTuv/9+/M3f/E3a+mqomILBYFr7iqxlRFIiIqJzw4i/fEdEROcOJiUiIvIMJiUiIvIMJiUiIvIMJiUiIvIMJiUasSoqKvDWW2+lOwwisoFJiYiIPINJiYiIPINJiUa8/v5+rF69GtOmTcO0adOwevVq9Pf3AwD279+PK6+8Es8++yz+9m//FtOmTcNLL72U5oiJMheTEo14zzzzDP70pz9h27Zt2L59Ow4cOID169fH2o8ePYrjx49j3759WL16Nerr6/HZZ5+lMWKizMWkRCPeK6+8grvuugtjx45FYWEh7rrrLmzfvj3WnpWVhbvuugvZ2dm46qqrkJeXh5aWljRGTJS5mJRoxOvs7ERJSUlsuqSkBJ2dnbHpgoICZGX99X6X5513Hvr6+lIaIxENYlKiEa+oqAjt7e2x6UgkEhsVmoi8hUmJRrzrr78ezzzzDLq7u9Hd3Y1169bhhhtuSHdYRGQgy3oWonPbwoULceLECcyePRsAMGvWLCxcuDDNURGREd5PiYiIPIOX74iIyDOYlIiIyDOYlIiIyDOYlIiIyDOYlIiIyDOYlIiIyDOYlIiIyDOYlIiIyDOYlIiIyDP+H4v2G2cP2+cnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "climatology.z.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "climatology.to_netcdf(f'{PREDDIR}climatology_{res}.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Climatology by week\n",
    "\n",
    "We can create amuch better climatology by taking the seasonal cycle into account. Here we will do this by creating a separate climatology for every week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "def create_weekly_climatology_forecast(ds_train, valid_time):\n",
    "    ds_train['week'] = ds_train['time.week']\n",
    "    weekly_averages = ds_train.groupby('week').mean('time')\n",
    "    valid_time['week'] = valid_time['time.week']\n",
    "    fc_list = []\n",
    "    for t in valid_time:\n",
    "        fc_list.append(weekly_averages.sel(week=t.week))\n",
    "    return xr.concat(fc_list, dim=valid_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "weekly_climatology = create_weekly_climatology_forecast(train_data, valid_data.time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (lat: 32, lon: 64, time: 17520)\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -87.19 -81.56 -75.94 -70.31 ... 75.94 81.56 87.19\n",
       "  * lon      (lon) float64 0.0 5.625 11.25 16.88 ... 337.5 343.1 348.8 354.4\n",
       "    week     (time) int64 52 52 52 52 52 52 52 52 52 52 ... 1 1 1 1 1 1 1 1 1 1\n",
       "  * time     (time) datetime64[ns] 2017-01-01 ... 2018-12-31T23:00:00\n",
       "Data variables:\n",
       "    z        (time, lat, lon) float32 dask.array<chunksize=(1, 32, 64), meta=np.ndarray>\n",
       "    t        (time, lat, lon) float32 dask.array<chunksize=(1, 32, 64), meta=np.ndarray>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weekly_climatology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "weekly_climatology.to_netcdf(f'{PREDDIR}weekly_climatology_{res}.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# The same for higher resolutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "for res in ['2.8125','1.40625']:\n",
    "    DATADIR = f'/media/rasp/Elements/weather-benchmark/{res}deg/'\n",
    "    # Load the entire dataset\n",
    "    z500 = xr.open_mfdataset(f'{DATADIR}geopotential/*.nc', combine='by_coords').z.sel(level=500)\n",
    "    t850 = xr.open_mfdataset(f'{DATADIR}temperature/*.nc', combine='by_coords').t.sel(level=850)\n",
    "    data = xr.merge([z500.drop('level'), t850.drop('level')])\n",
    "    # Load the validation subset of the data: 2017 and 2018\n",
    "    z500_valid = load_test_data(f'{DATADIR}geopotential', 'z')\n",
    "    t850_valid = load_test_data(f'{DATADIR}temperature', 't')\n",
    "    valid_data = xr.merge([z500_valid, t850_valid])\n",
    "    # Persistence forecast\n",
    "    persistence = []\n",
    "    for l in lead_times:\n",
    "        persistence.append(create_persistence_forecast(valid_data, int(l)))\n",
    "    persistence = xr.concat(persistence, dim=lead_times)\n",
    "    print(persistence)\n",
    "    persistence.to_netcdf(f'{PREDDIR}persistence_{res}.nc')\n",
    "    # Climatology\n",
    "    train_data = data.sel(time=slice(None, '2016'))\n",
    "    climatology = create_climatology_forecast(train_data)\n",
    "    print(climatology)\n",
    "    climatology.to_netcdf(f'{PREDDIR}climatology_{res}.nc')\n",
    "    # Weekly climatology\n",
    "    weekly_climatology = create_weekly_climatology_forecast(train_data, valid_data.time)\n",
    "    print(weekly_climatology)\n",
    "    weekly_climatology.to_netcdf(f'{PREDDIR}weekly_climatology_{res}.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# The End"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
